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Abstract 

The next-to-next-to-next-to-leading order (N 3 LO) Hamiltonian of potential non- 
relativistic QCD is derived. The complete matching of the Hamiltonian and the 
contribution from the ultrasoft dynamical gluons relevant for perturbative bound- 
state calculations is performed including one-, two-, and three-loop contributions. 
The threshold expansion is used to disentangle and match contributions of different 
scales in the effective-theory calculations. As a physical application, the heavy- 
quarkonium spectrum is obtained at N 3 LO for the case of vanishing QCD beta 
function. Our results set the stage for a full N 3 LO analysis of the heavy-quarkonium 
system. 

PACS numbers: 12.38.Aw, 12.38.Bx 



1 Introduction 



The theoretical study of nonrelativistic heavy-quark-antiquark systems |L| and its appli- 
cations to bottomonium || and top-antitop || physics rely entirely on first principles 
of QCD. These systems allow for a model-independent perturbative treatment. Non- 
perturbative effects @|| are well under control for the top-antitop system and, at least 
within the sum-rule approach, also for bottomonium. This makes heavy-quark-antiquark 
systems an ideal laboratory to determine fundamental parameters of QCD, such as the 
strong-coupling constant a s and the heavy-quark masses m q . The study of ti thresh- 
old production should even allow for a precision study of Higgs-boson-induced effects. 
Recently, essential progress has been made in the theoretical investigation of the nonrel- 
ativistic heavy-quark threshold dynamics based on the effective-theory approach 0J^| . In 
such a framework, one has two expansion parameters, a s and the relative velocity v of 
the heavy quarks. The corrections are classified by the total power of a s and v, i.e. N fe LO 
corrections contain terms of O (a l s v m ^, with l + m = k. This has the consequence that, in 
general, different loop orders, which are counted in powers of a s , contribute to the N fe LO 
result. 

Analytical results for the main parameters of the nonrelativistic heavy-quark-antiquark 
system are now available through the next-to- leading order (NLO) and the next-to-next- 
to- leading order (NNLO) [^,p|, p!0| , pT|p^ . |r^ . p^ ,|T5|] . They have been applied to bottomonium 
§[TT|rr|,[TB|rr^,|T3 and top-antitop [0|T§00,^g,gT|,0,g phenomenology. Some specific 
classes of the next-to-next-to-next-to-leading order (N 3 LO) corrections have also been in- 

(see Ref. 



vestigated [24,25,26,21 



28 



for a brief review). These corrections have turned 
out to be so sizeable that it appears to be indispensable to gain full control over this 
order, both with respect to phenomenological applications and in order to understand the 
structure and peculiarities of the nonrelativistic effective theory. Besides its phenomeno- 
logical importance, the heavy-quarkonium system is very interesting from the theoretical 
point of view because it possesses a highly sophisticated multiscale dynamics and its study 
demands the full power of the effective-theory approach. No qualitatively new theoretical 
effects are expected beyond N 3 LO, so that the N 3 LO analysis would bring us much closer 
to the full understanding of the perturbative nonrelativistic dynamics. 

In the present paper, we set the stage for the complete N 3 LO analysis of perturba- 
tive heavy quarkonium. In particular, we elaborate in detail the nonrelativistic effective 
Hamiltonian, which is the key object of the heavy-quarkonium theory, and its matching 
to the contribution associated with the emission and absorption of dynamical ultrasoft 
gluons, which are relevant for perturbative bound-state calculations [^4[^]. To this end, 
we employ the technique of Ref. [JJIJ based on the effective-theory approach gj^j31|] im- 



plemented with the threshold expansion [32 



This paper is organized as follows. In Section |2], we introduce the basic ingredients 
of the nonrelativistic effective-theory formalism and describe the main features of our 
approach. In Section [|, we recall lower-order results and analyze the one- and two-loop 
contributions to the N 3 LO Hamiltonian. In Section f|, we discuss the matching of the 
Hamiltonian to the ultrasoft contribution, which necessitates the inclusion of one-, two-, 
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and three-loop contributions. In Section |5], we convert our results into corrections to the 
heavy- quarkonium spectrum and illustrate their phenomenological relevance. Section |] 
contains our conclusions. In the Appendix, we present some details of our calculation of 
the l/m q corrections to the heavy-quark potential. 



2 Effective theory of nonrelativistic heavy quarks 

The nonrelativistic behavior of the heavy- quark- ant iquark pair is governed by a compli- 
cated multiscale dynamics. In the nonrelativistic regime, where the heavy-quark velocity 
v is of the order of the strong-coupling constant a s , the Coulomb effects are crucial and 
have to be taken into account to all orders in a s . This makes the use of the effective 
theory mandatory. The effective-theory approach allows us to separate the scales and to 
implement the expansion in v at the level of the Lagrangian. Let us recall that the dy- 
namics of a nonrelativistic quark-antiquark pair is characterized by four different regions 
and the corresponding modes 



(i) the hard region (the energy and three- momentum scale like m q ); 

(ii) the soft region (the energy and three- momentum scale like m q v); 

(hi) the potential region (the energy scales like m q v 2 , while the three-momentum scales 
like m q v); and 

(iv) the ultrasoft region (the energy and three-momentum scale like m q v 2 ). 

The ultrasoft region is only relevant for gluons, ghosts, and light quarks. Nonrelativis- 
tic QCD (NRQCD) ||[7| is obtained by integrating out the hard modes. Subsequently 
integrating out the soft modes and the potential gluons results in the effective theory of 
potential NRQCD (pNRQCD) which contains potential heavy quarks and ultrasoft 
gluons, ghosts, and light quarks as active particles. The effect of the modes that have been 
integrated out is two-fold: higher- dimensional operators appear in the effective Hamilto- 
nian, corresponding to an expansion in v, and the Wilson coefficients of the operators in 
the effective Hamiltonian acquire corrections, which are series in a s . 

The theory of pNRQCD is relevant for the description of the heavy-quarkonium sys- 
tem. Let us recall its basic ingredients. In pNRQCD, the (self)interactions between 
ultrasoft particles are described by the standard QCD Lagrangian. The interactions of 
the ultrasoft gluons with the heavy-quark-antiquark pair are ordered in v by the multipole 
expansion. For the N 3 LO analysis, only the leading-order (LO) emission and absorption 
of ultrasoft gluons have to be considered. They are described by the chromoelectric dipole 
interaction, which is of the form g s r ■ E, where g s is the QCD gauge coupling, r is the 
difference of coordinates of the quark and antiquark, and E = E a t a is the chromoelectric 
field, with t a being the generators of the colour gauge group SU(3). The propagation of 
the quark-antiquark pair in the colour- singlet (s) and colour-octet (o) states is described 
by the nonrelativistic Green function G s, ° of the Schrodinger equation, 

(H s >° - E) G s '°(r, r', E) = S(r - r'), (1) 
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where E is the energy of the quark-antiquark pair counted from the threshold 2m q and 
7i s, ° is the effective nonrelativistic Hamiltonian, 

US,o njS,o . 

it = T\q + • ■ • , 

n S c° = -— + VS>°(r), (2) 

Tflq 

with A r = d% and r = \r\. The ellipses stand for higher-order terms in a s and v. 
The Coulomb (C) potentials for the singlet and octet states are attractive and repulsive, 
respectively, and are given by 

VS(r) = -CjA 
r 

Voir) = " C F ) * (3) 

where = 3 and Cp = 4/3 are the eigenvalues of the quadratic Casimir operators 
of the adjoint and fundamental representations of the colour gauge group, respectively. 
Throughout this paper, we assume that a s = a s (fj) if no argument is specified. 

The LO approximation for the Green function is given by the Coulomb solution, 
which sums up terms singular at threshold and describes the leading binding effects. The 
corrections to the Coulomb Green function due to higher-order terms in the effective 
Hamiltonian can be found in Rayleigh-Schrodinger time- independent perturbation theory 
as in standard quantum mechanics. The Green functions have the following spectral 
representations: 

r s (r r > ^.f gWW i f d3fc gyjjW m 

where ip^ and are the wave functions of the quark-antiquark bound and continuum 
states, with principal quantum number n and relative three-momentum k, respectively, 
and the E + ie rule is implied. In Eqs. (|4]) and (|5p, the orbital and spin quantum numbers, 
I and m, respectively, are suppressed. Note that a discrete part of the spectrum (bound 
states) only exists for the singlet Green function. Through the emission or absorption of 
an ultrasoft gluon, the quark-antiquark pair changes its colour state, so that one switches 
from Eq. (|j) to Eq. (|5]) and vice versa. 

Let us now turn to the problem of perturbative calculations in the effective theory. 
Both NRQCD and pNRQCD have specific Feynman rules, which can be used for a sys- 
tematic perturbative expansion. However, this is complicated because the expansion of 
the Lagrangian corresponds to a particular subspace of the total phase space. Thus, in a 
perturbative calculation within the effective theory, one has to formally impose some re- 
strictions on the allowed values of the virtual momenta (see, e.g., Refs. |3^j3^] for examples 
of highly sophisticated calculations performed in this scheme). 
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Explicitly separating the phase space introduces additional scales to the problem, such 
as momentum cutoffs, and makes the approach considerably less transparent. A much 



more efficient and elegant method is based on the expansion by regions [[32,35], which is 
a systematic method to expand Feynman diagrams in any limit of momenta and masses. 
It consists of the following steps: 

(i) consider various regions of a loop four-momentum k and expand, in every region, 
the integrand in Taylor series with respect to the parameters that are considered to 
be small there; 

(ii) integrate the expanded integrand over the whole integration domain of the loop 
momenta; and 

(iii) put to zero any scaleless integral. 

In step (ii) , dimensional regularization, with d = 4 — 2e space-time dimensions, is used to 
handle the divergences. In the case of the threshold expansion in v, one has to deal with 
the four regions and their scaling rules listed above. 

In principle, the threshold expansion has to be applied to the Feynman diagrams of 
full QCD. However, after integrating out the hard modes, which corresponds to calcu- 
lating the hard-region contributions in the threshold expansion, it is possible to apply 
step (i) to the diagrams constructed from the NRQCD and pNRQCD Feynman rules |30| . 
Equivalently, the Lagrangian of the effective theory can be employed for a perturbative 
calculation without explicit restrictions on the virtual momenta if dimensional regular- 
ization is used and the formal expressions derived from the Feynman rules of the effective 
theory are understood in the sense of the threshold expansion. In this way, one arrives at a 
formulation of effective theory with two crucial virtues: the absence of additional regulator 
scales and the automatic matching of the contributions from different scales. The second 
property implies that the contributions of different modes, as computed in the effective 
theory, can be simply added up to get the full result. This automatic-matching property 
of effective-theory calculations in dimensional regularization was observed in Ref. and 
used for high-order calculations in the theory of QED bound states in Ref. []37| . We should 



emphasize, however, the crucial role of the threshold expansion in effective-theory per- 
turbative calculations because, in general, the naive use of the effective-theory Feynman 
rules and dimensional regularization leads to incorrect results. Another advantage of the 
effective-theory realization of the threshold expansion is that the individual contributions 
from the hard, soft/potential, and ultrasoft regions are manifestly gauge invariant. In- 
deed, the Lagrangians of NRQCD and pNRQCD are gauge invariant, and dimensional 
regularization preserves the gauge symmetry as well. Thus, the QCD calculation of the 
hard corrections and the NRQCD calculation of the soft and potential corrections to the 
on-shell amplitudes can be performed in the covariant gauge suitable for relativistic prob- 
lems, while the pNRQCD calculations can be done in the Coulomb gauge appropriate for 
nonrelativistic problems. In the next sections, we illustrate the power of the approach 
outlined above by an explicit analysis in pNRQCD at N 3 LO. 
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3 Nonrelativistic effective Hamiltonian 



Let us start this section with a general remark on the structure of higher-order corrections 
in pNRQCD. As already mentioned in the Introduction, the corrections are classified by 
the total power of a s , counting the number of loops, and v or, equivalently, l/m q . In 
particular, the N 3 LO Hamiltonian includes one-loop corrections of 0(a s v 2 ), two-loop cor- 
rections of 0(a 2 v), and three-loop corrections of 0(a 3 ). In NLO, the only source of 
corrections is the renormalization and running of the Coulomb potential. In NNLO, rela- 
tivistic corrections in v due to higher- dimensional operators start to contribute. In N 3 LO, 
retardation effects, which cannot be described by operators of instantaneous interaction, 
enter the game. They will be discussed in the next section. At this order, the nonrelativis- 
tic effective Hamiltonian becomes infrared (IR) sensitive. No qualitatively new theoretical 
effects are expected beyond N 3 LO. 

The general form of the Hamiltonian valid up to N 3 LO reads 



H = (2tt) 3 %) ( - -^V) + C c (a s )V c (\q\) + C 1/m {a s )V 1/m {\q\) + ^ ' h '' K 



m q 4m 3 



mi 



C s (a s ) + C P K) p2 + f' + C s (a s )S 2 + C x (a s )A(p, q) + C t (a s )T(q) 
Zq 



(6) 



where the following operators are involved 



Vc{\q\) 

Hp, q) 



2n 2 C F a s 



q 2 

,S-(pxq) 



m q \q\ 



Vi/ m (\q\) 
T(q) = (7 1 a 2 



(q ■ (Tx){q ■ tr 2 ) 



(7) 



Here, p and p' are the three-momenta of the incoming and outgoing quarks, respectively, 
q = p' — p is the three-momentum transfer, and cr 12 are the quark and antiquark spin 
operators. Note that the effective Hamiltonian is defined for on-shell quarks, with p 2 = 



p' 2 = m q E. The Wilson coefficients are power series in a s , 



CAa s 



E 

n=0 



Cn(m q , \q\,n), 



where the modified minimal-subtraction (MS) scheme for the renormalization of a s is 
implied. 

In the following, we discuss the nontrivial terms in Eq. (^) , which are sorted in terms of 
the inverse heavy-quark mass. The contribution from the hard-virtual-momentum region 
is analytic in v 2 and starts at 0{y 2 ). Thus, it does not affect Vcr(|<Jf|) and Vi/ m (|qr|). 
Corrections to the static Coulomb potential only arise from the soft contribution. Using 
renormalization group (RG) arguments, they can be rewritten as 



C c (a s )V G (\q\) 



47rC F a s (|q|) 



1 



<»«(|q|) 

47T 



ai + 



' a s (\q\y 
Air 



«2 
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An J \ q 
The RG logarithms can be recovered from Eq. (|9]) by recalling that (see, e.g., Ref. ||38|| ) 
a s (\q\) a 8 (fi) 



(9) 



7T 



7T 



1 + ! T* l+ ( 2 f) £ (« i + A 



+ 



'^(aO 



7T 



L^ 3 L 2 + ^o/9i^ + /9 2 ) + 



(10) 



where L 

aV 
aV 
aV 



= ln(/i 2 /g 2 ) and [[§ 
1 /2857 



64 



— Ci - ^Cft^n, - ™C A C F T Fni + 2CpTprii + ^-C A T F n* 



44 „ ^2 2 

¥ F fUi 



are the first three coefficients of the QCD beta function. Here, Tp = 1/2 is the index of 
the fundamental representation, and n\ is the number of light-quark flavors. The non-RG 
logarithmic term of O(o%) in Eq. (P) reflects the IR divergence of the static potential 
IfjQfl . The corresponding pole is canceled against the ultraviolet (UV) one of the ultrasoft 
contribution fl24|]4"l]l . For convenience, this pole is subtracted in Eq. (§) according to the 
MS prescription, so that the coefficient 03 is defined in the MS subtraction scheme both 
for UV and IR divergences. For consistency, the UV pole of the ultrasoft contribution has 
to be subtracted in the same way. Throughout the calculation, we use the same procedure 
to render the contributions from the various regions finite. The actual cancellation of the 
spurious divergences appearing in the process of expanding by regions is reflected in the 
(jl independence of the final result. 

In the literature pi,H|4T]1, the coefficient in front of the IR logarithm in the 0(a^) 
static potential is given as C^/(247r), which differs from C a /(Stt) in Eq. @. This is a 
consequence of the consistent use of dimensional regularization in our analysis based on 
the threshold expansion. The difference is due to the fact that we perform all three loop 
integrals in d dimensions, not just the one that is IR divergent. The logarithmic terms not 
associated with IR-divergent integrals are unphysical and are exactly canceled by similar 
terms from the three-loop ultrasoft-potential-potential contribution, in which only the 
ultrasoft integral is U V divergent (see Section |]) , while the physical logarithmic integral 
between soft and ultrasoft scales results in lna s corrections to the spectrum. The calcu- 
lation of the coefficients a« can be performed in the static limit m q — > 00 of NRQCD. Due 
to the exponentiation of the static potential [j42 |, these coefficients only receive contribu- 
tions form the maximum non-Abelian parts. In the language of the threshold expansion, 
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the selection of these parts effectively separates the contribution of the soft region. The 
Abelian colour factor Cp indicates the presence of the Coulomb singularity and implies 
that at least one loop momentum is potential. All such contributions are just iterations 
of the lower-order potential and are taken into account in the perturbative solution of the 
Schrodinger equation ([!]) around the Coulomb approximation. 
The one-loop coefficient, 

ai = yOA - —T F m, (12) 

has been known for a long time |42],43|, while the two-loop coefficient has only recently 
been found {|4],f|5]. In our previous communication ||30|| , we confirmed the result of 
Ref. 



15 



a 2 



4343 



162 



+ 47T 2 



- + -C(3) 
4 3 SV ' 



r 2 



1798 



C A T F ni 



f " 16C(3) 



C,,T r i)i ( yT F n ( 



(13) 



where ( is Riemann's zeta function, with value £(3) = 1.202057. . .. At present, only Pade 
estimates of the tree-loop MS coefficient are available, namely ||6 



03 

4 3 




if ni = 3, 
if ni = 4, 
if ni = 5. 



(14) 



Although the Pade estimates are in reasonable agreement with the exact results where the 
latter are available, the reliability of Eq. ([14] ) is not guaranteed, and it is very desirable 
to exactly evaluate the coefficient 0,3. However, in Section |5|, we show that even a 100% 
uncertainty in a 3 would not result in a significant error in the N 3 LO corrections to the 
spectrum for the states with small principal quantum number. 

The l/m^-suppressed terms of Eq. (|j) receive contributions from the soft and potential 
regions and, applying RG techniques, can be written in the following form: 



Ci/ m (a s )Vi /m (|q|) = — 



m q \q\ 



71 



b 2 -Uc 2 A + 2C A C F )ln^ 



+ 



The one-loop coefficient b\ reads [47|48|] 



-C A + 



c f 



(15) 



(16) 



Our result for the two- loop coefficient b 2 is listed in Ref. [3C]. A more detailed description 
of the calculation is presented in Section ^.2. The coefficient of the two-loop IR logarithm 
in Eq. ( |i~5l) can be extracted from the UV divergence of the ultrasoft contribution [24 . 



However, similarly to the case of the IR divergence of the static potential, one has to take 
into account additional logarithmic terms resulting from the consistent use of dimensional 
regular ization. 
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The RG analysis results in the following representation of the l/m q part of the effective 
Hamiltonian: 



•6 , a s(v) 



(]P _|_ aS ^ \ J? 



71 



a s (n)C 5 (a s ) = a B (\q\) < d + 

rv 

a s (n)C p (a s ) = a s (\q\) 
a s (n)Ci(a s ) = a s (\q\) 
a s {n)Ci{a s ) = a s (m q ) 



d{ + -(C A -2C F )ln^ 



+ 



d\--C A \n^\ 



4 + ^d[ + 



7T 



'/, 



a . a s( m q) da 



TV 



q 

i — s,X,t 



i = 5,s, 



(17) 



where the contributions from the annihilation channel are marked by the superscript 
a. The normalization scale of a s in the one-loop scattering terms is not fixed because 
they receive contributions from both the soft and hard regions. According to the RG, 



the scale of a s should be chosen to be \Xh 



for the hard contribution and ji s 



a s m q for the soft one. The IR one-loop logarithms which match the UV behavior of the 



ultrasoft contribution |24| are written out explicitly in Eq. (0). The calculation of the 
one-loop coefficients d\ within the threshold expansion is discussed in Section |3|.l. Some 
of these coefficients also contain logarithms of the form In (m 2 q /q 2 ^ originating from the 
logarithmic integration between the soft and hard scales. 

The purely relativistic tree-level 0{y 2 ) corrections are given, up to a colour factor, by 
the standard Breit Hamiltonian and read 



d'i 



d s 

5, a 



(7 s 
"o 



rfn =6, 



dl 



1 

3' 



0. 



d s ' a = 0. 



The QED effective Hamiltonian for n; light fermions is obtained from the above ex- 
pressions by setting C A = and Cp = T F = 1. Note that the QED Breit Hamiltonian 
has a nonvanishing one-photon annihilation coefficient, 



d, 



0,QED 



1. 



(19) 



which is absent in the case of colour-singlet quarkonium due to colour conservation. 
In Refs. 



36 



dealing with QED bound-state calculations in the Coulomb gauge 
one finds Wilson coefficients different from the Abelian parts of those listed in Eqs 
and fli~8|) . Using them would lead to 



dn 



1. 



(20) 



These two sets of coefficients are equivalent, and the difference is related to the use of off- 
shell operators in the Hamiltonian. This problem is discussed in more detail in Section 13.2. 
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3.1 One- loop operators 

The 0(a s v 2 ) operators, contributing to N 3 LO at one loop, have attracted some attention 
and an essential part of the results can be found in the literature [^, ^J49|J50|J5T|J5^J53|J54 



Here, we present a consistent derivation of these corrections within the threshold-expansion 
framework. The 0(a s v 2 ) operators receive contributions from the hard and soft/potential 
regions. The contribution from the hard region requires a fully relativistic treatment. A 
part of it is directly related to the Wilson coefficients Cp, cp, eg, and d 2 parameterizing the 
Fermi, Darwin, spin-orbital, and heavy-quark vacuum-polarization terms in the NRQCD 
Lagrangian, respectively ||52|| . The residual part arises from the on-shell scattering and 
annihilation box diagrams at threshold |xj. 

The calculation of the soft contribution can be performed in NRQCD. In the effective- 
theory language, we study the reduction from NRQCD to pNRQCD and compute the 
effect of the soft modes being integrated out. Apart from the standard LO terms of the 
NRQCD Lagrangian, we need the l/m g -suppressed terms originating from the covariant- 
derivative operator D 2 / (2m q ) acting on the quark and antiquark fields, the Fermi, Darwin, 
and spin-orbital terms. Note that the covariant- derivative operator includes the quark 
kinetic-energy term k 2 /(2m q ). According to the threshold expansion, it should be treated 
as a perturbation if k is soft or kept in the nonrelativistic quark propagator, 

S W = i — vnk ^ - » (21) 

k — k /{2m q ) + is 

if k is potential. The potential region is connected with the contribution of the pole of 
the propagator of Eq. ( pT|) to the integral over k . In this connection, we can make an 
interesting observation. Let us consider two-particle-irreducible diagrams that are free 
of Coulomb singularities, so that the ko contour is not pinched between the poles of 
the quark and antiquark propagators. In this case, one can close the contour of the ko 
integral keeping the poles either inside or outside. The contribution from the potential 
region is obviously different in these cases, although the result for the integral is the 
same. This means that separating the contributions from the soft and potential regions 
for diagrams without Coulomb singularity is useless, and the propagator of Eq. (pi] ) can 
thus be safely expanded in l/m q in both regions. In fact, the expansion of the pole 
contribution yields familiar generalized functions of k , namely 5^(k ). This observation 
dramatically simplifies the calculation, which can be performed in a covariant form after 
substituting k = vo ■ k, where vq = (1, 0). 

By contrast, the two-particle- reducible diagrams including the product of the quark 
and antiquark propagators, 

(22) 



k — k I (2m q ) +ie k + k / (2m q ) — is 

where k is the two-particle-reducible loop four-momentum, suffer from a Coulomb singu- 
larity. In this case, after expanding the quark propagator, one obtains ill-defined pinched 
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products like 



1 1 

(fco + ie) m (fc - ie) 7 



(23) 



Thus, separating the soft and potential regions is unavoidable. In the soft region, the pole 
contributions of the quark and antiquark propagators have to be excluded, and the above 
product should be defined to be its principal value, 



1 



(fc + is) 



+ 



1 



(k 



is 



\m+n 



(24) 



which again allows for a covariant treatment. In the potential region, the quark and 
antiquark propagator poles produce contributions of the form 



VK 



k 2 



IE 



6 h 



k 2 



2m n 



+ 5[ko + 



2m„ 



(25) 



where the 1/v Coulomb singularity shows up explicitly. After integration over /c , Eq 
yields the nonrelativistic Green function of the free Schrodinger equation. A contribution 
of this type can always be related to iterations of the operators of the effective Hamiltonian 
in time-independent perturbation theory. Therefore, it should not be considered as a 
correction to the Hamiltonian. However, there is one subtle question here, namely as to 
whether the operators that vanish for on-shell quarks should be included in the effective 
Hamiltonian or not. This does not affect the potential, but it matters for the l/m q 

operator discussed in Section [|2. 

One has to be careful with the definition of commutators of the Dirac/Pauli matrices 
within dimensional regularization. Since poles in e are only present in the individual 
contributions from the hard and soft regions and drop out in the sum, one can explicitly 
retain the commutators during the analysis of the these regions and replace them by the 
four /three-dimensional expressions in the final result ]37[. Otherwise one has to use the 
same prescription for the evaluation of the commutators of the Dirac/Pauli matrices in 
the hard and soft regions. Throughout the calculation, we use the four /three-dimensional 
antisymmetric e tensor for the definition of the commutators as was done in Ref. 
(This differs from the prescription of Ref. |53|j .) 



54 



The one-loop calculation poses no technical problems, and our results for the Wilson 
coefficients read 




4 + -(c A 



+ 



i.n4 

3 mi 



d\ - -C A In —z 



C A + -^-T F ni 



C f — — Tp 
15 
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d\ 

d\ 
d\ 



d, L 



■S.a 



4 7 u 2 



mi 









) ° A - 3^ 


+ 


[( 




ft 





4 + 21n 



3 



6 m„ 



C A + 4C F 



+ 



^-4 + 41n2-i27r)T F , 
(2 - 2 In 2 + ztt)T f , 



14 7 
In — 

27 6 q 2 , 



6 o 2 



13 



C A - Y Tpni 



108 6 



In 



Ca T F m 

27 



(26) 



where the contributions from the hard (h) and soft (s) regions are explicitly separated. 
The first two equations of Eq. (^) are written in a way appropriate for Eq. (|H]). As in 
NNLO, the one-photon-annihilation channel provides an additional contribution to the 
QED Wilson coefficient, namely 



d\ 



s,a 
QED 



2 In 2 + ITT - ( - 



2 III. 2 ■ i« ) y. 



(27) 



The QED part of Eq. (^) is in agreement with Ref . |JT| . The non-Abelian part of Eq. 
agrees with Ref. ||54|| , including the nonlogarithmic Ca term in the coefficient d\, which 



differs from the one of Ref. 48 



3.2 Two-loop operators 

The O {a 2 s v ) part or the effective Hamiltonian is given by the two-loop corrections to 
the l/m q potential. These corrections are solely generated by the covariant-derivative 
operator in the NRQCD Hamiltonian. The calculation of the two-loop l/m q potential is 
simplified by the absence of the hard contribution, but, in turn, it is complicated by the 
presence of off-shell operators. 

To introduce the problem, let us start with the one- loop l/m q potential. For illus- 
trative purposes, we use the Feynman gauge, where the Coulomb and transverse gluons 
do not mix. Nonzero contributions come from the planar and crossed box diagrams and 
a diagram with one three-gluon vertex. The last two diagrams do not contain Coulomb 
singularities. According to the procedure described in the previous section, the quark and 
antiquark propagators should be expanded in l/m q . The result reads 

-(C A C F -C F )^. (28) 

The soft contribution to the planar box diagram vanishes because the l/m q terms from 
the expansion of the quark and antiquark propagators cancel. The potential contribution 
corresponds to the second iteration of the operators generated by the exchange of one 
potential gluon. The operators that are defined for on-shell quarks enter the effective 
Hamiltonian, and their iteration is taken into account when Eq. ([I]) is solved. However, 
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the potential-gluon exchange also generates operators proportional to the energy transfer 
Qo = (p' 2 ~P 2 )/(^ m q)^ which vanish for on-shell quarks. Such operators do not enter the 
effective Hamiltonian, which is defined on-shell, and their iterations should be considered 
as corrections to the effective Hamiltonian. The factor go cancels the denominator of 
the free nonrelativistic Green function and effectively makes the diagram two-particle 
irreducible. The Coulomb singularity brings a factor of m q , so that for the calculation of 
the l/m q corrections we need the iteration of the 1/m 2 off-shell operator and the leading 
Coulomb potential. In the case under consideration, the relevant off-shell operator is 



TcC F CX s fp' — p 2 



m 2 \ q 2 



(29) 



Explicit evaluation of the corresponding potential contribution yields 



£3^2 (30) 

2m q \q\ 



and the coefficients of Eqs. ( Pq) and (|30|) sum up to Eq. (|lq) . In fact, by using the 
Coulomb equation of motion, it is straightforward to check that the matrix elements of 
Eqs. (|29f) and fl30D between Coulomb states are the same. 



Note that, in QED calculations performed in the Coulomb gauge p6| , |49"| , the off-shell 
tree-level operator analogous to Eq. (^) naturally appears. If one includes this operator in 
the effective Hamiltonian, the Wilson coefficient d s Q is given by Eq. ( P0"D and the coefficient 
bi is purely non-Abelian. 

The equivalence of the two formulations is obvious from the above analysis. The use of 
off-shell operators is advantageous in QED because it allows one to reduce the number of 
loops by means of the Coulomb equation of motion, as may be seen by comparing Eqs. ( |2"9"| ) 
and (pop. However, here we use the on-shell formulation and the general covariant gauge, 
which is more suitable for multiloop QCD calculations. 




Figure 1: Examples of two-particle-irreducible two-loop diagrams. The standard quark-gluon 
vertex represents the leading Coulomb interaction. The black circles correspond to the three 
types of 0(l/m q ) terms generated by the quark covariant-derivative term. 



The structure of the expansion remains intact at two loops, and our final result reads 
b 2 = - + | In 2) C\ + - ~ In 2) C A C F + ^C A T Fni - \c F T F n x . (31) 
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There is no fully Abelian Cp contribution in 62, as is known from the QED analysis |36|J^9|] . 
In the calculation of the two-loop two-particle-irreducible diagrams, which completely 
determine the maximum non- Abelian C\Cf and CaCfTfUi structures of the result, we 
used an expanded form of Eq. (|2l|). Typical diagrams contributing to the C\Cp part are 
depicted in Fig. [L|. The analysis of the relevant two-particle-reducible diagrams, shown 
in Fig. D, is conceptually similar to the one-loop case described above. Let us discuss 
it in more detail. The reducible loop momentum can be either soft or potential. If it 
is soft, the quark and antiquark propagators can be expanded, and the only nonzero 
contribution corresponds to the situation where the single gluon is the Coulomb one and 
the l/m q term is kept in the expansion of the one-loop block, B. If the reducible loop 
momentum is potential, then one only has to take into account contributions from the 
off-shell operator. There are two possibilities: (i) the off-shell operator is generated by 
the single-gluon exchange, and the block B stands for the one-loop corrections to the 
Coulomb potential; or (ii) the off-shell operator comes from the l/m^ part of the one- 
loop block B, and the single gluon is the Coulomb one. The analysis of the potential 
contribution is rather straightforward from the technical point of view. The calculation 
of the two-particle-irreducible diagrams and the soft parts of the two-particle-reducible 
diagrams is more involved. Some details are presented in the Appendix. 




Figure 2: Example of a two-particle-reducible two- loop diagram. B stands a general one-loop 
two-particle-irreducible subgraph. The threshold expansion of this diagram is discussed in the 
text. 



We performed a number of nontrivial checks for our analysis, (i) We worked in the 
general covariant gauge and verified that the gauge parameter cancels in our final result, 
(ii) The two-loop expression from which Eq. ( |3l"D is obtained contains both UV and IR 
divergences. The UV ones were removed in Eq. (^Tj) by the renormalization of a s in the 
one-loop result of Eq. (|T5|). On the other hand, the IR divergences were canceled by 
the UV ones of the ultrasoft contribution (see Section f|) leaving a finite result for the 
spectrum. The RG logarithms proportional to (3q and the IR logarithms are in agreement 
with Eq. (]15D. (iii) To test our program, we also recalculated the two-loop corrections 
to the static heavy-quark-antiquark potential and found agreement with Ref. P5] . Note 
that, with our prescription for the calculation of the soft and potential contributions, we 
explicitly obtained zero for the partially Abelian corrections to the static potential. 



14 



4 Ultrasoft contribution 



For the N 3 LO Hamiltonian, only the leading retardation effects are needed. They arise 
from the chromoelectric dipole interaction of the heavy quarkonium with a virtual ultrasoft 
gluon, as depicted in Fig. [| In the analysis of the ultrasoft contribution, we proceed along 
the lines of the original analysis ^4|. As has been mentioned in Section [^, there is freedom 
in the choice of gauge. We work in the Coulomb gauge, which is especially appropriate 
for N 3 LO calculations in pNRQCD because the Coulomb gluon does not propagate and 
the dynamical gluon is transverse. 




Figure 3: Feynman diagram giving rise to the ultrasoft contribution at N 3 LO. The shaded and 
light double lines stand for the singlet and octet Green functions, respectively. The loopy line 
represents the ultrasoft-gluon propagator in the Coulomb gauge, and the black circles correspond 
to the chromoelectric dipole interaction. 



The analytical expression for the corrections to the singlet Coulomb wave function can 
be obtained by using the pNRQCD Feynman rules of Section |2|. It reads 

SG\x,y,E) = -^ r „^y: (y) n ^), (32) 

where 

r d 3 k ( k 2 \ 

J mn {E) = -C F g]{^) J -^{ rj ) mk { ri ) kn I» [E-—J, (33) 

The sum/integral in Eq. (^) goes over the whole spectrum, and m and n stand for 
the complete set of quantum numbers characterizing the discrete/continuum part of the 
spectrum. The scale fi ns ~ a 2 s m q of g s in Eq. (^) reflects the ultrasoft-momentum flow 
through the gluon propagator. The matrix element (r)k n is taken between the singlet 
Coulomb wave function of quantum number n and the octet Coulomb wave function of 
three-momentum k. Performing in Eq. Q3~4] ) the integration over l Q , we recover the ex- 
pression of time-independent nonrelativistic perturbation theory. The remaining integral 



with 
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over I is UV divergent. Subtracting the UV pole according to the MS prescription, we 
obtain 



I ij (E-k 2 /m q 



5 ij 
6tt 2 



E-—) 3 I In 



c 



k z /m q - E 



+ ln T ^ T + ^-ln2 



(35) 



where, for convenience, the Coulomb energy, 



E 



c 



C F a 2 m q 
4n 2 ' 



(36) 



with n — 1, has been introduced into the arguments of the logarithms. The fc-dependent 
logarithmic term in Eq. (|35|) represents a pure retardation effect and cannot be interpreted 
in terms of some instantaneous interaction. It receives contributions from Coulomb- 
gluon-exchange diagrams of all orders and leads to a QCD analogue of the familiar Bethe 
logarithms in the corrections to the spectrum, to be discussed in the next section. On the 
other hand, making use of the completeness relation, 



d 3 fc 

(2vr) 3 



Je-—) 1 = (r(E-H° c ) l r) 
\ m q J \ fr 



(37) 



the remaining part of Eq. (fffi), excluding the fc-dependent logarithmic term, is reduced 
to an instantaneous interaction of the form 



r{E-n° c fr 



^^-(C\ + 2C A C F 



+ Ca — ^ (Ay-, - \) + reducible part, 
mi I r J / 



(38) 



where the reducible part includes terms with the operator E — 7i s c acting directly on 
a wave function. By using Eqs. (|32|) , (|35|), (^7|), and (|38D , one arrives at the following 
representation of the corrections to the Green function: 

5G s (r, r', E) = - J d 3 r"G s c (r, r", E)H us (r")G s c (r" , r', E) + contact terms, (39) 

where, in momentum space, 



C F a s 



1, V 2 



(E?Y 



In 2 



C^ + 4(C 2 A + 2C A C F 



m q \q\ 



a x 



-8(C A -2C F )^ + 1QC A ^ 



a s p 2 +p' 21 



m- 



2q 2 



(40) 



The first term in Eq. ( |3"9"| ) imitates the corrections to the Green function due to the 
N 3 LO term of Eq. (|40"D in the effective Hamiltonian. The contact terms correspond to the 
reducible part of Eq. (|38l). By the equation of motion ([]]), one or both Coulomb Green 
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functions in Eq. (^) are converted into S functions. The corresponding contribution 
cannot be imitated by a term of the effective Hamiltonian. It does not affect the energy 
levels, but it leads to corrections to the wave functions. A part of such corrections, namely 
the on-shell renormalization of the heavy-quarkonium wave function at the origin, was 
computed in Ref. p^fl . In the present paper, we refrain from discussing this type of 
corrections. 





(5»- 




Figure 4: Examples of two- and three-loop diagrams encoded in the diagram of Fig. (3|, which 
require additional matching to bring Eq. (^) in agreement with the threshold expansion. The 
dashed and loopy lines represent the potential (Coulomb) and ultrasoft (transverse) gluon prop- 
agators in the Coulomb gauge, respectively. The black circles correspond to the interaction 
generated by the quark covariant-derivative term. 



Although dimensional regularization was used in deriving Eq. ([40]), the latter is not 
consistent with the threshold expansion because the three-dimensional expression for the 
Coulomb Green function was used instead of the (d — l)-dimensional one, which is not 
available. Thus, Eq. (fh]) derived in Ref. ]24| requires some additional matching. For 
this purpose, we separate the divergent contributions from the diagram of Fig. ^ and 
compute them according to the threshold expansion. Eq. (fHf ) implies that only the 
one-, two-, and three-loop contributions encoded in the diagram of Fig. |3| are divergent. 
They include one divergent ultrasoft loop integration and zero, one, or two convergent 
potential loop integrations. The three-dimensional form of the Coulomb Green function 
used in the evaluation above implies that the integration over the potential momenta 
is performed in three dimensions, while, for the corresponding contributions obtained 
within the threshold expansion, they are done in d — 1 dimensions. Thus, the matching 
term is given by the difference of the diagrams computed in d — 1 dimensions and the 
same diagrams with three-dimensional integrals over the potential momenta in the limit 
e — > 0. Only two- and three- loop diagrams have to be considered. Examples of such 
diagrams are presented in Fig. |j. The calculation is simplified by the fact that, for the 
matching, we only need the pole part of the ultrasoft integral, which factorizes. For 
the two-particle-reducible diagrams of the type shown in Fig. |2], one has to take into 
account the off-shell operators generated by the two-particle-irreducible block B with one 
ultrasoft loop in a way similar to the case of the l/m q potential in Section |3|.2. In addition 
to the 0(a s v 2 ) operator proportional to Eq. ( p9|) with an extra factor of a s , the one-loop 
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ultrasoft exchange generates an off-shell operator proportional to 



irC F a s P 2 + P' 2 - 2m q E 
m 2 q q 2 



(41) 



which, in three dimensions after adding one extra Coulomb gluon, results in the same 
l/m q term as the operator of Eq. (|29|) . The matching terms are found to be 



Incidentally, the three-loop coefficient does not have a constant term. The \i dependence 
of the ultrasoft contribution TC ns + 5H ns exactly cancels the \i dependence of the N 3 LO 
Hamiltonian given in Eqs. @, (|15|) , and (|17| ) ensuring the cancellation of IR and UV 
poles. 

5 Heavy- quarkonium spectrum 

Let us now apply the results of the previous sections to the analysis of the heavy- 
quarkonium spectrum. We restrict the analysis to the perturbative corrections, neglecting 
issues like nonperturbative contributions in the case of bottomonium and finite-width ef- 
fects in the case of the top-antitop system. This is justified because the problem of 
large perturbative corrections seems to be crucial for the heavy-quarkonium theory. Fur- 
thermore, we only consider the zero-orbital-momentum states, with I = 0, which are of 
primary phenomenological interest. 

5.1 Perturbative cx^m q heavy-quarkonium spectrum 

The O (a 3 ) corrections to the energy levels arise from several sources: 

(i) matrix elements of the N 3 LO operators of the effective Hamiltonian between Cou- 
lomb wave functions; 

(ii) higher iterations of the NLO and NNLO operators of the effective Hamiltonian in 
time-independent perturbation theory; 

(iii) matrix elements of the N 3 LO instantaneous operators generated by the emission 
and absorption of ultrasoft gluons; and 

(iv) retarded ultrasoft contribution. 




(42) 
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Parts (i) and (ii) include corrections due to the running of a s in the lower-order operators 
of the effective Hamiltonian proportional to /3j, with i = 0,1,2. The logarithmic part of 
these corrections can be taken into account by choosing the relevant soft normalization 
scale of a s in the NNLO result for the spectrum to be p s ~ a s m q . However, there are 
also nonlogarithmic corrections proportional to the QCD beta function. We postpone 
the calculation of these corrections to a future publication. Here, we focus our attention 
on the conceptually interesting non-RG corrections. In the absence of the running of 
a s , the calculation of part (ii) is reduced to a redefinition of a s in the leading Coulomb 
approximation. The matrix elements relevant for parts (i) and (iii) are conveniently 
evaluated in coordinate space. All necessary formulae, including Fourier transforms, can 



be found in Ref. [48|. Part (iv) corresponds to the fc-dependent term of Eq. (^). The 
ultrasoft corrections to the n-th energy level are given by J nn (E n ), and its retarded part 
can be written as 

9r ,3 n/ 3 

' (43) 



3tt 

where we have introduced the QCD Bethe logarithms |24 



^ :\(r)kn\ 2 ( E° - ^Vln - E ) 2 - . (44) 



n CWMJ (2tt)^ V n mj EC-k 2 /m q 

The latter can be reduced to one-parameter integrals of elementary functions H24jl.fl For 
the reader's convenience, we list the relevant formulae here. They read 

L E n = / dvY n E (v)X*(v), (45) 
Jo 



where 



yE/ \ _ 2 6 p^z/(z/ 2 + 1) exp[4z/arccot(z//p n )] n 2 v 2 



n 2 (i> 2 + p 2 l ) 3 [exp(27ru) — 1] v 2 + p 2 n 

X 1 (u)=p 1 + 2, 

Yf v is 2 (2p 2 + 9p 2 + 8)-p 2 (p 2 + 4) 
[y 2 + pi) 



u\8pl + 60p 2 + 123p 3 + 66) - 2u 2 p 2 (6p 2 + 41p 3 + 54) + 3p|(p 3 + 6) 

XM ~ W r T P W ' (46) 



Pn = n(^--l) = ^. (47) 



with 

Pn — n, 

\ZLjf / o 

The expressions for X n with n > 3 are usually irrelevant for practical applications. For 
n = 1, 2, 3, we obtain the following numerical values: 

Lf = -81.5379, Lf = -37.6710, Lf = -22.4818. (48) 



1 There are two misprints in Ref. P4fl : in Eq. (A. 3), n should be in the numerator; in Eq. (A. 5), arctan 
should be replaced by arccot . 
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Putting everything together and writing 

E n = E? + 5E^ + 5E® + 6E® + ■■■ 



(49) 



we obtain our final result for the N 3 LO corrections to the heavy- quarkonium energy levels 
in the approximation of putting (3{a s ) = 0: 



5E^ 





E C n 




P(a s )=0 







a x a 2 + a 3 
32tt 2 



+ 



CaCf 
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4 16n 2 



Cl 
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r 3 



n 



139 7 41 
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V108 12n 6 U as V v 

79 7 8 7 

7^-^r+Q ln2 + - Inn + ^ x (n + 1)) + dL Qs - 
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n 

5(5 + 1)' 



d2 



15 

49C A C F r F ^ 



2 In 2 + (1 - In 2)5(5 + 1) 



C F T F 



n 



36n 



1 ^5(5 + 1) 



9 18n 27 



C F T F m 2 



71 



n 



(50) 



where 5 is the spin quantum number, L as = — \n(C F a s ), ^i(z) = dT(z)/dz, T(z) is Eu- 
ler's gamma function, and je — 0.577216 ... is Euler's constant. The terms proportional 
to ai correspond to iterations of lower-order operators. We have not included in Eq. (|50|) 
the imaginary part corresponding to the partial width of the decay of the 5 = state to 
two gluons, 
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The logarithmic part of Eq. (|50| ) , 

„3 rl 



6E^ 



E 



c 



a: 



s 

7T I 6 



C\ + — C F C\ 4 
in 



C F T F a b s m q 
2n~ 3 



41 7 2 

6n 6n 3n 2 



(51) 



C 2 F C A + -C% \ In 
n 



0:.s 



(52) 

is known from previous analyses p5| , |2"o| . For the corrections to the n = 1,2,3 energy 
levels, we find numerically 



5E{ 



(3) 



a 



c 



E\ 



Q3 

d27T 2 



177.716 - 11.611m + 0.274 nf - 0.004 nf 
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6E^ 



SEP 



60.500 In — + f— 18.853 + 1.312 n { + 6.222 In — ) S(S + 1) 

a. \ a. 



71 



E. 



c 



- 33.389 In 



-^r + 102.917 - 8.034 m + 0.274 n] - 0.004 n 3 

1 



7T 



E. 



c 



«3 



-9.603 + 0.658 n t + 3.111 In —)S(S + 1) 

OL s ) 

75.919 - 6.690 n x + 0.274 n] - 0.004 nf 



_32tt 2 

- 23.957 In — + f— 6.425 + 0.439 n { + 2.074 In — ) S(S + 1) 

a s \ a s 



(53) 



5.2 Numerical estimates and phenomenological examples 

To illustrate the phenomenological relevance of our results, let us consider two important 
physical examples: the resonance in top-antitop threshold production by e + e~ annihi- 
lation via a virtual photon and the lowest T resonance. We neglect nonperturbative 
contributions and finite-width effects, so that the masses of the resonances are deter- 
mined by the perturbative expressions with principal quantum number n — 1 and spin 
quantum number £7 = 1. The complete NLO and NNLO corrections may be found in 
Refs. [3, P^ , P!8| , and the N 3 LO ones for (3{a s ) = are given in Eq. fl5H|). To take into 
account the N 3 LO RG logarithms, we normalize a s in NLO and NNLO at the soft scale 
fA s = CFa s (fi s )m q . The setting of the normalization scale in the C(a 3 ) corrections is a 
more subtle problem. In this order, the hard and ultrasoft regions start to contribute. 
This results in RG logarithms with corresponding scales. Furthermore, the contributions 
from different regions are not separately finite, and the operators of the effective Hamil- 
tonian acquire anomalous dimensions, which result in non-RG logarithms [see Eqs. @, 
(Hjl)) and (I7D1- Thus, starting with the next order, the RG and non-RG logarithms mix. 
The correct treatment of the logarithmic corrections is possible within the effective-theory 
RG approach [PD| , ^| , ^ , ^7| , ^j5*^|j6l]|j6l| . For simplicity, we ignore these sophistications for 



the time being and employ the soft normalization point for the whole of the 0{a s ) cor- 



rections. As an estimate of the nonlogarithmic corrections proportional to 

1 



we use the 



0k term, which is currently only known for n 
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1 7l Z 7T 4 

"8 + 16 + 1440 



+ C(3)-yC(3) + ^C(5) 



(54) 



with £(5) = 1.036928 . . ., and is expected to dominate the corrections of this type. 

For the top-antitop system, we use m t = 176 GeV and a s (fi s ) = 0.14 to find the 
perturbative expansion of the binding energy to be 



E x = -1.53 GeV x [1 + 0.448 + 0.322 + (0.006 + 0.011| a3 + 0.073^) + 



(55) 



where the N 3 LO contributions due to the and 0q terms are given separately. Thus, 
the resonance peak position is decreased by approximately 140 MeV in comparison to 
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the NNLO result. Although the 0(a z s ) corrections are still important, the series shows 
reasonable convergence, which makes us optimistic about an accurate determination of 
rrit and a s from this observable. 

As for the T(15) resonance, we use nib = 4.8 GeV and a s (/x s ) = 0.31 to find the 
perturbative expansion of the binding energy to be 

Ex = -205 MeV x [1 + 1.11 + 1.88 + (0.49 + 0.19| a3 + 1.02^) + ■■■]. (56) 

This implies that the value of m& extracted from the T(1S) resonance is increased by 
approximately 170 MeV in comparison to the NNLO result. Although there is no further 
growth of the perturbative corrections, the 0(a z s ) corrections seem to be too large to 
expect a reliable prediction from the N 3 LO result, and some optimization, e.g. by mass 
and/or coupling-constant redefinition, is necessary to improve the convergence of the 
series. 

In the above estimates we used the Pade results of Eq. (ITJj) for the coefficient a 3 . The 
accuracy of the Pade approximation is difficult to estimate, and a significant deviation 
from the exact result does not seem impossible. However, from Eqs. fl55l) and (|56|) we 
observe that the corresponding contribution only provides about 10% of the total C(af) 
corrections, and even a 100% variation of 03 merely results in a 10% variation of the 0(oP s ) 
corrections. This is not crucial for the top-antitop system, but it could be essential in the 
bottomonium case, where the magnitude of the 0{oP s ) contribution is very sizeable. The 
analytical evaluation of the coefficient is thus quite important. 

The origin of the large NLO and NNLO corrections in Eqs. ( [55]) and (|56|) is usually 
attributed to the IR-renormalon contribution. This contribution is absent in the IR- 
safe "short-distance" masses, and the perturbative series for such masses are expected 
to exhibit faster convergence (see, e.g., Refs. ||i5| , |ll||i7l1 ). We observe, however, that, in 
the top-antitop case, the perturbative series numerically converges even in the pole-mass 
scheme. On the other hand, in the bottomonium case, the 0(a^) corrections remain 
sizeable even for /3(a s ) = 0, i.e. the naive subtraction of the renormalon contribution 
through a mass redefinition does not completely solve the problem of large perturbative 
corrections. 

Our next comment concerns the corrections logarithmic in a s . Using the effective- 
theory RG equations, it is possible to sum up the logarithmic corrections to the energy 
levels to all orders. The presence of several correlated scales renders the problem very 
interesting and nontrivial from the conceptual point of view. Now there are two contra- 
dictory results on the resummation of the a™ +4 ln n a s terms in the series for the 
energy levels, the first of which is given by Eq. fl5"2"|) . The result of Ref. [^0| is obtained 



within pNRQCD, while an alternative formulation of effective theory, velocity NRQCD 
(vNRQCD), was used in Ref. |p9|. Since there can be only one correct result, this issue 
has to be clarified. In any case, it is interesting to check the accuracy of the logarithmic 
approximation. Numerically, the logarithmic series is dominated by its first term, which 
provides about 80% in both the bottomonium and top-antitop cases |5^|J6T1|J52| . From 



Eq. (|53"D, we observe that the full result has approximately the same magnitude, but the 



opposite sign compared to the logarithmic contribution. We thus conclude that, while the 
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logarithmic contributions in some order can give us a hint at the order of magnitude of 
the full contribution of that order, the practical relevance of the high-order resummation 
is questionable. This is not unexpected because the resummation parameter a" ln m a s is 
neither large for a s ~ nor for a s ~ 1. 



6 Conclusion 

In this paper, we took a crucial step towards the N 3 LO analysis of the heavy-quark thresh- 
old dynamics. We used the effective theory of pNRQCD and the threshold expansion for 
a detailed analysis of the nonrelativistic Hamiltonian in this order. Explicit expressions 
for the N 3 LO Hamiltonian in one and two loops were given. We also presented the full 
matching of the Hamiltonian to the contribution from the ultrasoft gluons, which enter 
the bound-state dynamics in this order. The matching calculation includes one-, two-, 
and three-loop operators. To complete our analysis, the three-loop MS coefficient 03 of 
the corrections to the static potential, for which only Pade estimates are available, has to 
be computed. 

With the full expression for the Hamiltonian and the ultrasoft contribution at hand, 
it is straightforward to complete the N 3 LO analysis of the heavy-quarkonium spectrum. 
In this paper, we derived the heavy-quarkonium spectrum in this order, neglecting the 
nonlogarithmic terms proportional to the QCD beta function. For the latter, only the 0q 
term for the ground-state energy is known so far. 

Collecting all available contributions, we found the N 3 LO corrections to be sizable 
for the top-antitop system, where, however, the perturbative series for the resonance- 
peak energy exhibits a tendency to converge even in the pole-mass scheme. In the case 
of the T(IS') resonance, the corrections are so sizeable, that some kind of optimization 
of the perturbative expansion is needed, e.g. by removing the pole mass in favour of a 
renormalon-free short-distance mass. However, we should emphasize that, in the bot- 
tomonium case, the N 3 LO corrections remain sizeable even for (3(a s ) = 0, and the bad 
behaviour of the perturbative series cannot be solely explained by the renormalon contri- 
bution. 

In order to render the analysis more accurate, the remaining nonlogarithmic terms 
proportional to the QCD beta function have to be evaluated along with the 03 coefficient. 
However, we do not expect a qualitative change of our result. 

The result of this paper also provides a starting point for the calculation of the N 3 LO 
single-logarithmic a 3 In a s corrections to the heavy-quarkonium production and annihila- 
tion rates, similarly to a number of QED bound-state problems p3| , |64| , |65| , |66| . While the 
analysis of the nonlogarithmic terms in this order requires the calculation of the three- 
loop hard renormalization of the relevant production and annihilation amplitudes, which 
still is a challenging theoretical problem, the logarithmic terms are universal and essen- 
tially determined by the effective Hamiltonian and the ultrasoft contribution presented 



in this paper. Along with the known double-logarithmic a 3 In a s terms [26|, the single- 



logarithmic contribution would constitute an essential part of the N 3 LO corrections to 
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the heavy- quarkonium production and annihilation rates. The calculation of the single- 
logarithmic terms would also lead to further progress in the resummation of the logarithms 
in v via the nonrelativistic effective-theory RG |58U6T| , since it determines an anomalous 
dimension necessary for the NNLO logarithmic analysis of heavy- quarkonium production 
and annihilation. Although it is, in general, dangerous to rely on the logarithmic ap- 
proximation, as we observed in the case of the spectrum, the situation could be different 
for the cross section normalization, where the N 3 LO double-logarithmic contribution is 
known to be sizeable and the resummation of the logarithmic terms could stabilize the 
behavior of perturbation theory []5g] . 
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Appendix 



To evaluate the two-loop contribution to the heavy-quark potential, in particular the l/m q 
corrections, one needs analytical results, at least as Laurent expansions in e up to some 
order (typically, e° and e 1 ), for the following family of two-loop Feynman integrals: 

2 N r r d d kd d l 
J(oi,...,a 8 ;g ;e) 



(h?)*i(l2)aa[(k- q y]<H[(l 
1 



q) 2 } a *[(k - /) 2 ] ffl 5 



x 



(57) 

(v ■ k) a z (v o ■ l) a7 [v ■ (k - l)] as ' 

where the four- vector vq is defined below Eq. (pl~D, k and I are loop four-momenta, and 
+ie is omitted in all the denominators. 

As in Ref. [[44] , fi5| , |B"7|l , we use a reduction procedure that expresses any integral of the 
form of Eq. (^7|) in terms of some master integrals. To this end, we employ the following 
identities, which can be obtained by means of the method of integration by parts f6"8fl , 
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as well as the trivial identity 6 — 7 = 8 . Here, the standard notation for raising and 
lowering operators has been used, e.g. 



1 _ 3 + J(ai, 



,a 8 j 



J(ai — 1, a 2 , a 3 + 1, 



a 8 1 



(59) 



We developed a reduction procedure very similar to the one of Ref. |57| . In our problem, 
however, we need a larger class of integrals that arise in the calculation of the l/m q 
corrections in the general covariant gauge. The main difference between our reduction 



procedure and the one of Ref. |67j is that we stop the reduction if we arrive at integrals 
expressed in terms of gamma functions for finite e. There are two subclasses of the 
integrals of Eq. (|57|) that are only evaluated as expansions in e up to some order. The 
first of them was described in Ref. fUjl , namely 



J(ai,...,o B ;g ; e) 
In particular, we have 



d d kd d l 
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J(0, a 2 , a 3 , 0, a 5 , a 6 , 0, a 8 ) = 7(a 5 , a 3 , a 2 , a 8 , a 6 ), 
J(ai, 0, 0, a 4 , 05, a6, «7, 0) = /(a-i, 04, 05, a6, 07). 

The master integrals for this subclass are 
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(62) 



We need also integrals of the type of Eq. (|60| ) with a numerator that can be chosen to be 
q-k or q-l. The reduction of these integrals is quite similar and also results in two master 
integrals. 

The second subclass of the integrals of Eq. (^) that are not expressed in terms of 
gamma functions for general e consists of integrals with 05 = = a-j = 0. Using 
Feynman parameters, the general integral of this kind can be represented in terms of the 
following Mellin-Barnes representation: 

(_l)a 2 «6-l ( l7{ d/2\ 2 
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(63) 
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where a = Yh=i a %- We performed a reduction to integrals with 

and evaluated them by means of Eq. (|63|) for the required values of a%. In particular, we 
found 



J(l, 1, 1, 1, 0, 0, 0, 1; q 2 ; e) = \ ^ /7 J e { -An 2 In 2[1 + e(4 + In 2)] + — e + 0(e 2 ) 



: _ 9 2)i/2 



(64) 
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